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We propose a versatile variational method to investigate the spatio-temporal dynamics of one- 
dimensional magnetically-trapped Bose-condensed gases. To this end we employ a g-Gaussian trial 
wave-function that interpolates between the low- and the high-density limit of the ground state 
of a Bose-condensed gas. Our main result consists of reducing the Gross-Pitaevskii equation, a 
nonlinear partial differential equation describing the T = dynamics of the condensate, to a set of 
only three equations: two coupled nonlinear ordinary differential equations describing the phase and 
the curvature of the wave-function and a separate algebraic equation yielding the generalized width. 
Our equations recover those of the usual Gaussian variational approach (in the low-density regime) , 
and the hydrodynamic equations that describe the high-density regime. Finally, we show a detailed 
comparison between the numerical results of our equations and those of the original Gross-Pitaevskii 
equation. 



I. INTRODUCTION 

Due to the unprecedented experimental maneuverabil- 
ity of Bose-condensed gases, their intrinsically nonlinear 
dynamics became an appealing research topic engrossing 
scientists from fields as diverse as numerical analysis and 
condensed matter physics P, [3, [1| ■ The T = dynamics 
of these gases is accurately described by the so-called 
Gross-Pitaevskii equation (GPE), a cubic Schrodinger 
equation brought forward in the early sixties by two sem- 
inal papers of Gross and Pitaevskii 4, 5]. 

While there is good quantitative agreement between 
the full GPE and current experimental results on a wide 
range of topics, including: dark soliton dynamics in- 
side maganetic traps 0, Ql, bright soliton interactions 
Q, non-equilibrium oscillations in binary BECs and 
Faraday waves in periodically driven BECs there 
are, however, rather few models giving analytical in- 
sight. Naturally, reducing the partial-differential, there- 
fore infinitely- dimensional, GPE to the level of a few 
ordinary differential equations (ODEs) is a nontrivial 
task. Previous works include the so-called Gaussian vari- 
ational model (which describes the low-density regime) 
[O, d; El) [i3l and the hydrodynamic analysis of high- 
density condensates [l^, [ll] . All these works addressed 
the dynamics of a Bose-condensed gas close to its ground 
state. For a comprehensive account, from both an ex- 
perimental and theoretical viewpoint, of excited states 
such as solitons, vortices and alike we refer to the re- 
cent review 17] and references therein. In this paper 
we present an efficient universal variational model able 



to describe both the high- and the low-density dynamics 
close to the ground state {i.e., without entering the realm 
of soliton- like solutions). Our model stems from the q- 
Gaussian trial wave-functions that were used to describe 
the ground-state properties of a trapped Bose-condensed 
gases [18|, ll9|, |20j. Though restricted to ground-state 
properties, the work of Fetter [ll] is particularly rele- 
vant since he is the first one to point out the universality 
of this ansatz. 

The rest of the paper is structured as follows: In Sec.HIl 
we present the variational method using the g-Gaussian 
trial wave-function. In Sec. Ill Al we reduce the varia- 
tional equations to the level of three coupled ODEs and 
in Sec. Ill Bi we simplify these equations to the level of two 
coupled nonlinear ordinary differential equations describ- 
ing the phase and the curvature of the wave-function, 
and a separate algebraic equation yielding the generalized 
width. In Sees. IIIII and IIVI we show that our equations 
recover those of the variational Gaussian model and the 
standard hydrodynamic formulation of the GPE. Section 
IVl presents our numerical results for a particular case of 
a GPE with a periodically forced nonlinearity. Finally, 
Sec. IVIl is alloted to our concluding remarks. 



II. VARIATIONAL MODEL 

A. Coupled equations 

The equation describing the T = dynamics of a 
trapped three-dimensional Bose-condensed gas is given, 
in adimensional units, by 
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where the time-dependent nonHnearity is given by U{t) = 
47ra(t), a(t) is the time-dependent two-body scattering 
length, and we have taken for simphcity h — m = 1. 

In this work we consider a quasi-one-dimensional 
Bose-Einstein condensate (BEC) confined in a strongly 
anisotropic magnetic trap. In such a setup the dynam- 
ics of the condensate is well-known to follow a one- 
dimensional version (with rescaled nonlinearity coeffi- 
cient) of the three-dimensional GPE ([1]) [21]. Now that 
the system under scrutiny is one-dimensional in nature, 
we can consider the usual magnetic, parabolic, trapping 
potential of the form 



Vix) 



(2) 



Let us now use the following trial wave-function (from 
now on called the g-Gaussian) to describe the state of 
the condensate close to its ground state 



^(x,<)-/[q(t)]Viv(^l-ii 
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The g-Gaussian ansatz (jS]) has three, free, time- 
dependent parameters {w,(3,q} corresponding, respec- 
tively, to the width (which is strictly positive at all 
times), chirp, and q parameter, where the q parameter 
captures the regime of the BEC cloud (see below). To 
determine / [q{t)] one imposes the normalization condi- 
tion 



\tlj{x,t)\^ dx = N, 



(4) 



D 



where N represents the total number of atoms in the 
condensate and 



D 



(5) 



is the domain supporting the g-Gaussian ansatz. Com- 
puting the norm ([¥]) one finds 
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where B(-, •) is the usual Euler beta function (23 |. 

It is crucial to note the two important limits of the 
g-Gaussian ansatz: in the q(t) — > 1 limit our ansatz re- 
covers the usual Gaussian ansatz (valid for low-density 
condensates), as 



lim 



1 - 
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while the q{t) ^ — 1 limit recovers the usual parabolic 
density profile of the Thomas-Fermi region (valid for 
high-density condensates), as 



lim 
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and 
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The novelty of this ansatz is that in addition to having 
the usual variational parameters 'w{t) (width of the con- 
densate) and fi{t) (chirp), it also accounts for the possible 
changes of curvature through the q{t) variable. As the 
above limits clearly show q{t) is the crucial ingredient for 
being able to recover analytically both the low- and the 
high-density regime of a BEC. Naturally, the g-Gaussian 
is not the only function that is able to interpolate be- 
tween these two limits of a BEC [l^. The so-called 5'„ 
function 



Sn{x) = exp ( -2_] — 
\ fc=l / 



(11) 



becomes a Gaussian for n — 1 and turns into a parabolic 
profile for n oo (and \x\ < 1), as 



Soc{x) 



°° ^2fe 



exp 



fe=i 



exp (ln(l-a;2)) = 1 



The advantage of using the g-Gaussian lies however in the 
fact that this function is amenable to analytical compu- 
tation (and gives rise to very simple equations describing 
the dynamics of the condensate), while other functions 
(like the Sn{x)) are not. 

Our variational method follows the traditional recipe: 
we introduce the g-Gaussian trial wave-function in the 
BEC Lagrangian defined as 



h{t) = / C{x,t)dx, 



(12) 



where D is given by ([5]) and the Lagrangian density is 
given by 

(13) 

The dynamics of the BEC, restricted to our ansatz, is 
then determined through the Euler-Lagrange equations 



d 5L _ aL 
dt dy dy ' 
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where y = {w,/3, g}. These equations will be the main 
result of our paper and, as we will confirm in the fol- 
lowing sections, they are able to accurately describe the 
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are the four terms of the BEC Lagrangian. The only 
apparently cumbersome term is L2 which involves 



A[qit)] 
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(20) 



A close numerical inspection of this function shows that 
it is, however, almost linear, as can be easily seen from 
Fig.ffl 

For practical purposes one can use the two linear ap- 
proximations: 



FIG. 1; (Color online) The solid black line shows A [q(f)] as 
given by Eq. (|20|) . while the blue dashed line and the red 
dash-dotted line show A [q{t)] as given by the approximate 
formulas (ISTJ and (l22l). 



dynamics of the condensate in both the high- and the 
low-density limit. 

A couple of straightforward (though tedious) integra- 
tions show that 



where 



m 

N 



U{t) 

L3(t) 
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A[q{t)]^ 
for the low-density case and 
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A[q{t)] = 
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for the high-density one. These linear approximations 
come from the series expansion of A [q{t)] around q(t) — 1 
and —1. It is worth mentioning that in our numerics (see 
SeclYl we use the full expression ((20|) . 

Let us notice that by setting q{t) to —1 the first term 
in L4 diverges, a result which is perfectly consistent with 
the Thomas-Fermi theory, where one neglects the kinetic 
term, which is numerically small in the bulk of the con- 
densate for large number of atoms, due to its divergence 
at the border of the cloud. As will be shown shortly, our 
variational method works flawlessly in the high-density 
limit where q(t) will be close but always larger than — 1, 
therefore the divergence will never be an issue. As far as 
we know this is the first variational method that is able 
to describe dynamically the Thomas- Fermi limit. 

The Euler-Lagrange equations for {w, /3, q} are 



qit)~5 A[q{t)]NUit) , 2n^^^^^l + ml + mi ^ ^ ^33^ 



4ii;(t)3(l + g(t)) V2wit)^ " 7 - 3q{t) 

Ht) (^m) - 7^^) = 2^(^) (24) 
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I 

where {■)q stands for derivative with respect to q. Adding an extra exponential term like exp^icjj) in the 
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FIG. 2: The line shows F [q(t)] as given by Eq. To illus- 

trate the singularity at q(t) = — 1 we have chosen a log- plot. 
Please notice that outside the vicinity of q{t) = —1 F [q{t)] is 
numerically well behaved. Notice as well that F [q{t)] has an 
infinite slope at q{t) = ±1. 



trial wave-function, i.e., extending the set of variables to 
{N, (f>, w, P, q}, gives rise to N = 0, the equation assuring 
the conservation of the norm, and (f> = 0, but leaves the 
previous equations unchanged. We have dropped this 
term for simplicity. Naturally, one can imagine various 
other corrections to the present ansatz. However, even 
simple phase corrections such as exp (i7„(t)a;"), where 
n is an even integer strictly bigger than two, yield an 
analytically intractable Lagrange function. 

While Eqs. (P5)) -(|25 p are relatively simple to write 
down, they are in fact differential- algebraic equations and 
(if left in their current form) require specialized numeri- 
cal treatment (see, for instance, the monograph of Haircr 
et al. HI). 

To understand the nature of this dynamical system let 
us remember that variational methods are usually build 
around pairs of canonically conjugated variables., which 
give rise to coupled ODEs. Taking, for example, the sim- 
ple case of g = 1 one easily notices that p = /2 (the 
generalized momentum) and q = P (the generalized coor- 
dinate) are the two canonically conjugate variables. Our 
ansatz, however, does not posses this intrinsic symmetry, 
as q has no canonical conjugate, which finally gives rise 
to the algebraic constraint. Since the Lagrange function 
becomes analytically intractable at the slightest change 
in the ansatz, it is unfruitful to try working with pairs of 
conjugate variables. This latter option requires the use 
of a direct method [1^ , which has the same lack of trans- 
parency as the numerical solution of the original GPE. 
Instead, in order to gain analytical insight, one can easily 
decouple the previous equations, a task which is detailed 
in the next subsection. 



B. Decoupled equations 

In order to simplify the numerical solution of equations 
Eq. (ESll-dSni) one multiplies Eq. 1^ by 3w{t)/2{7-3q{t)) 
and then subtracts it from Eq. The ensuing equa- 

tion is purely algebraic. Its solution with respect to 'w{t) 
is 



w{t) = 



NU{t) ' 



(26) 



where 
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■3qit)) + A[q{t)]} 



(27) 



Equation (|26p can now be used in the other two equa- 
tions. We have now decoupled the equation of the width, 
w{t), from those of q{t) and /3(i). A close numerical 
inspection of F [q{t)] shows that it has a singularity at 
q(t) = — 1 (which, as mentioned before, follows physically 
from the divergence of the kinetic energy at the border of 
the cloud) , and is (numerically) well behaved everywhere 
else as one can easily infer from Fig. ^ 

The original variational equations ((23|) ~ ((25)) can now 
be recast as 



m 



'U{t) 



F,[q{t)] 3 ■ 

■ F [q{t)] 3q{t) - 7 



(28) 



for q{t) and 

m 

for /3{t), where 
G[q{t)] 



= -2/32-^+G[g(i)] N^Uit)\ 



(29) 



(5 - q{t) + 2V2A [q{t)] F [q{t)] (1 + g(t))) 
x{7-3q{t))F~* [q{t)]{l + q{t))-'/16 

{3A [qit)] + 2A, [qit)] (7 - 3qm' (7 - 3qf 
x{3A[qit)]-A[qit)] iqit)-5){l + qm 
x8{l + qit)fiq{t)-9)-\qit)-l)-^/81. 



Equations (pS)) and (|29p constitute the main theoretical 
result of our work. 

Before going into the high- and low-density limits of 
these two equations let us make the following comment: 
the dynamics of q{t) depends only indirectly on the num- 
ber of particles. It is therefore extremely instructive to 
analyze the denominator in Eq. 



Em] 



'F[q{t)] 3q{t)-7' 



(30) 



see Fig. 3, and notice that it diverges in the vicinities of 
q{t) = ±1 and is well behaved everywhere else. Let us 
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FIG. 3: (Color online) The line shows E [q{t)] The divergence 
around q{t) — ±1 indicates that the low- and the high-density 
regimes have a very weak surface dynamics, i.e., q{t) « 0. The 
red dot (of coordinates {qc,E[qc]) = (0.27605..., -7.76722...)) 
and the accompanying dashed red line (of equation q = qc = 
0.27605) indicate the extremum value of E [q{t)] at qc. 



notice that the extremum value of E [q{t)] corresponds 



to Qc « 0.27605.... The divergence of E[q{t)] originates 
from the infinite slope of F [q{t)] at q{t) — ±1 and in- 
dicates that q{t) w or, in other words, that in both 
the low- and the high-density regime the dynamics of a 
BEC affects only superficially the curvature of the wave 
function. This shows in turn why the use of a Gaussian 
ansatz in the Thomas-Fermi regime gives rise to quanti- 
tatively good results (see, for instance, Ref. The 
main point we want to stress is that these two physically 
different regimes have in common a very weak surface 
dynamics, therefore almost any ansatz with an intrinsic 
breathing mode provides decent quantitative results. 

Finally, let us notice that the Jacobi matrix of the dy- 
namical system made out of Eqs. and (com- 
puted with the ground state values of q and /3 and, 
of course, constant U) has complex conjugate roots 
throughout the whole density spectrum. 

III. HIGH-DENSITY CONDENSATES 

To see how the previous equations recover the usual 
hydrodynamic formulation let us consider for A[q(i)] the 
approximation given by Eq. (|22p . Within this approxi- 
mation we have that 
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(32) 



As F[q[t)\ reduces to 
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(33) 



and we are investigating the dynamics of the system in a 
vicinity of q{f) = —1, Eqs. (PT|) and ([5^ take the simpler 
form 
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(34) 



and 



(35) 



Equation ((35|) along with the Eq. (f24|) . where we now 
ignore the time derivative of q{f), i.e., 



w{t) = 2w{t)P{t), 



(36) 



are the so-called hydrodynamic equations. To see this 
more clearly let us work in the variables a(t) and a(t) 
defined as 



Pit) 
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and 



a{t)=- 



3iV 
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so that equations (|3ip and ([5^ read 

a{t) + 3a{t)a{t) 
a{t) + a{tf + Q"^ + 2U{t)a{t) 



(37) 



(38) 



(39) 
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These are exactly the hydrodynamic equations put for- 
ward by Dalfovo et al. |16| . 



IV. LOW-DENSITY CONDENSATES 

A. Gaussian ansatz 

In the low-density limit one can set q{t) to 1 and ig- 
nore its derivative with respect to time, in which case 
Eqs. ([23ll -([25 |l can be recast as 



therefore 
w{t) 



w{t) = 4w{t)(3{tf + 2w{t)P{t), 
1 NU{t) 



2TTw{ty 



+ n^w{t) = 0. 



(41) 



(42) 



This is the well-known equation that describe the dy- 
namics of the width of the condensate in the low-density 
limit. 



NU{t) n^w{t) 



2w(tf 2^/2^w(^)2 2 

+2w{t)l3{tf +w{t)l3{t) = 

2w{t)P{t) = w{t). (40) 

From the second equation in one has that 



B. g-Gaussian ansatz 

One way to refine the standard equation stemming 
from the Gaussian ansatz is to start from Eqs. (pS)) and 
and use for A [q{t)] the high-density approximation 
provided by Eq. (PT|). The ensuing equations are 



and 
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As F[q{t)] reduces to 
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and we are investigating the dynamics of the system in a 
vicinity of q{t) — 1, Eqs. and (|^^ take the simpler 

form 



C/(i) 



and 



NU{t) 



2 2w{ty u{q)V2^w{t)^'^ 



where 



u{q) = 



128 



142 - 81q{t) + 27q{tf 



(46) 



(47) 



(48) 



Replacing in Eq. ([^5)1 q{t) with 1 in all factors except 
the one yielding a zero term, one has that 



NU{t)w{t) 
4a/2^ 



(44) 



(49) 



Finally, inserting q{t) in Eq. one has that 

w{t) = 2w{t)(3(t). (50) 

Except for the u{q) factor Eqs. (|47p and (|50p are iden- 
tical to Eqs. (gOl)- Of course, the two sets of equations 
coincide in the N = limit, but otherwise the ones stem- 
ming from the g-Gaussian ansatz have the clear advan- 
tage of accounting for the change in the curvature of the 
wave- function. 

V. NUMERICAL RESULTS 

In the previous two sections we have shown that the 
equations stemming from the g-Gaussian ansatz recover 
those of the usual Gaussian approach (in the low-density 
regime) and the standard high-density hydrodynamic 
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FIG. 4: Dynamics of the BEG density profile under the mod- 
ulation of the nonlinearity for different atom numbers A*' (as 
indicated on each plot). Parameter values correspond to: 
n = f , (7o = f , e = O.f and uj = f .6. 



equations. We now turn to a detailed numerical com- 
parison between our variational equations and the orig- 
inal GPE. To this end, we consider a magnetic trap of 
the form ([2]) and a modulated nonlinear term U{t) = 
Uq{1 + tsm.{iot)). The modulation of the nonlinear term 
can be achieved by modulation of the transverse width 
(see Ref. [l^ and references therein) or by, the so-called, 
Feshbach resonance technique [1^, [l^ . 

For illustration purposes we choose a magnetic trap 
with SI = 1 and a nonlinearity modulation with Uq — 1, 
e — 0.1 and uj = 1.6. Let us notice that the natural 
frequency of the condensate is equal to Qcon, where ujn 
depends weakly on N and is numerically close to 2, there- 
fore choosing uj — 1.6 assures a series of beats (or, tech- 
nically, parametric resonances) between the frequency of 
the driving field and the natural frequency of the conden- 
sate as it can be observed in the space-time plots of the 
density for different values of N in Fig. [5] (for a more de- 
tailed analysis of mode-locking and energy transfer due to 
resonances see Refs. [28] and [29] and references therein). 

The advantage of our approach is that close to reso- 
nance a system responds strongest to excitations, there- 
fore this is the ideal setup to investigate the limits of our 
equations. In order to probe the different regimes of the 
condensate we choose five different cases corresponding 
to A'' = 0.1, 1, 10, 10^, 10'^. These cases encompass the 
bulk part of the density spectrum, the corresponding val- 
ues of q (restricted to only three significant digits) being 
equal to 0.990, 0.906, 0.431, -0.236, and -0.710, respec- 
tively. In Figs. [5] and [S] we depict, for various values of 
N, the steady state solution u{x) to the GPE alongside 
the best linear square fits using the g-Gaussian, the tra- 
ditional Gaussian (see Fig. [5]) and the parabolic Thomas- 
Fermi approximation (see Fig. [S]). The steady state u{x) 
of Eq. (ID) is obtained by taking 



ip{x, t) = u{x)e 



-ifit 
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where /i is the so-called chemical potential. The ensuing 
steady state equation is then solved using a fixed point 



FIG. 5: (Golor online) Comparison between the GPE steady 
state solution, the g-Gaussian fit and the Gaussian fit. The 
top four panels show the steady state solution to the GPE 
(thick black dashed line), its fit using the g-Gaussian (thick 
black solid line, almost indistinguishable from the dashed 
line due to the very good approximation) and the traditional 
Gaussian (thin red solid line) for different atom numbers A''. 
The bottom four panels show the respective differences be- 
tween the steady state and the g-Gaussian (thick black line) 
and the steady state and the traditional Gaussian (thin red 
line). 



iteration algorithm (Newton method) by optimizing /i to 
give the desired total mass. In Figs.[5]and[6]we also depict 
(see lower set of panels) bM(x) = ugpe(2;) — Mfit(2;), i.e., 
the difference between the actual steady state of the GPE 
and the corresponding fits. As it can be clearly seen from 
both figures, the g-Gaussian provides an extremely good 
approximation of the wave-function in hoth the low- and 
the high-density regimes. In fact, as it can be evidence 
from Fig. [51 even in the low-energy regime, where it is 
well known that the traditional Gaussian ansatz is a good 
approximation, the g-Gaussian is superior. Incidentally 
the same is true when analyzing the high-density regime 
where, as it can be evidence from Fig. [51 the g-Gaussian 
provides a better approximation than the Thomas-Fermi 
profile. 

To quantify the accuracy of our variational treatment 
of the GPE we will follow two distinct veins: i.) we 
fit the spatial component of the time-dependent den- 
sity profile provided by a GPE solver with a g-Gaussian 
profile and compare the ensuing time-series for [9(0] j 
w(t), and q(t) with the corresponding ones obtained from 
Eqs. (|23)) - (|25|) : also, we compute the full width at half 
maximum (FWHM), wfwhmIOi ^ function of time 
both for the density profile provided by the GPE solver 
and that obtained from Eqs. (|^5)) - ([^5)) : ii.) using two dif- 
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FIG. 6: (Color online) Comparison between the GPE steady 
state solution, the g-Gaussian fit and the parabolic Thomas- 
Fermi profile fit. The panels show the same information as in 
Fig. [5] by replacing the traditional Gaussian by the parabolic 
Thomas-Formi profile. 
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ferent fits based on a Gaussian and a Thomas-Fermi den- 
sity profile, and the g-Gaussian profile that yields from 
Eqs. ((23)) - ([25]) . we calculate the time-dependent mean 
square deviation (MSD) of these profiles from that ob- 
tained from the original GPE. 

Figs. Em and El depict the p [qit)], w(t), and q{t) 
time-series for, respectively, the low-density regime {N — 
0.1, 1), the high-density regime {N = 10^, 10^), and the 
intermediate density regime {N = 10). Notice the very 
good agreement for [q{t)] and w{t) for all the cases 
under consideration. 

In the low density limit, when the condensate cloud is 
close to a Gaussian profile, one can observe from Fig. [7] 
that the dynamics of the parameter q for the full GPE 
numerics and the reduced ODE system are in good agree- 
ment. However, in the high- and intermediate density 
limits, as it can be evidenced from Figs. [8] and [9l there 
seems to be a (somewhat counterintuitive) poor agree- 
ment in the q time-series despite the fact that [q{t)] 
and w(t) are still in very good agreement. This apparent 
shortcoming should not be seen as a fault of our varia- 
tional treatment, but rather as an indication of the fitting 
sensitivity of the g-component of our density profile. In 
fact, we have observed that even the addition of small 
quantities of noise (< 5%) to the profile to be fitted in- 
duce sizeable changes in the resulting q value (results 
not shown here). To understand this more clearly one 
should look into how the fitting procedure works: when 
minimizing the mean-square deviation between the GPE 
profile and the g-Gaussian one encounters the derivative 
of the g-Gaussian with respect to q which has a log-type 
of singularity at the border of the cloud. This singularity 



FIG. 7: (Color online) Low-density regime. Comparison be- 
tween the full GPE dynamics and the reduced ODE dynamics 
using the g-Gaussian approach. The top (bottom) series of 
panels corresponds to atom numbers A'' = 0.1 (A'' — 1). For 
each atom number, the three panels corresponds to, from top 
to bottom, the time-series of a) the peak density f^[q{t)], 
b) width w(t), and c) the g-Gaussian parameter q{t). Thick 
blue lines corresponds to the full GPE dynamics while thin 
red lines to the reduced ODE dynamics using the g-Gaussian 
variational approach. Note that for f^[q{t)] and w{t) the time- 
series obtained from the full GPE dynamics and our ODE 
reduction are practically indistinguishable. 

makes it very hard to make a one-to-one comparison of 
q(t). One way to avoid this problem is to compute q(t) 
(via Eq. ([26]) ) from the w{t). This latter approach (see 
Fig. [T0|) shows good agreement between the two time- 
series. 

Another way to test the quality of the g-Gaussian 
ansatz is to make a direct comparison between the 
WFWBMit) time-series (see Fig. [TT|). The t«FWHM(i) of 
the GPE density profile is computed numerically, while 
that of the g-Gaussian is given by 



2\/-2 + 2i/2(i+9(t)) 

WFWHuit) = , w{t). (51) 

^/-l + q{t) 

Please notice the trivial limits q{t) 1 and q{t) ^ — 1, 
in which cases we get the well-known results wfwhm (t) — 
2Vlog 2w{t) and Wfwhm(0 = V^w{t). 

In order to further substantiate the quality of the q- 
Gaussian approach we depict in Fig. [T^lthe mean square 
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FIG. 8: (Color online) High-density regime. All panels depict 
same information as in Fig. [7] but in this case for A'' — 10^ 
and iV = 10^ 
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FIG. 9: (Color online) Intermediate density regime. All pan- 
els depict same information as in Fig. [7| but in this case for 
N = 10. 



deviation between the different fits (note the logarith- 
mic scale in the vertical axis). In the low-density regime 
the g-Gaussian function performs some four orders of 
magnitude better than the Thomas- Fermi fit. This re- 
sult is rather obvious since the Thomas-Fermi limit is 
not the correct limit to use in this case. However, when 
comparing the g-Gaussian (blue solid line) with the tra- 
ditional Gaussian (red dashed line) one can still notice 
more than one order of magnitude better agreement for 
the g-Gaussian trial function. On the other extreme of 
the regime (i.e., the high-density regime), the g-Gaussian 
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FIG. 10: (Color online) Comparison of the g-Gaussian param- 
eter q{t) between full GPE dynamics (thick blue line) and the 
reduced ODE dynamics by fitting q{t) from w{t) (thin red 
line). Each panel, from top to bottom, corresponds to: a) 
N = lQ-\ b) A- = 10", c) A = 10^ d) N = 10^ and e) 
N = 10''. 
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FIG. 11: (Color online) Full width at half maximum 
'ii'FWHM(i). Each panel, from top to bottom, corresponds to: 
a) A = 10~\ b) A = 10°, c) A^ = 10^ d) A^ = 10^ and 
e) TV = 10''. Thick blue lines corresponds to the full GPE 
dynamics while thin red lines to the reduced ODE dynamics 
using the g-Gaussian variational approach. 



performs about two orders of magnitude better than the 
traditional Gaussian. This is again obvious since the tra- 
ditional Gaussian is not the right limit to use in this 
regime. On the other hand, it is interesting to note that 
the g-Gaussian has a similar performance to the Thomas- 
Fermi approximation in this high-density regime. Finally, 
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FIG. 12: (Color online) Mean square deviation (MSD) be- 
tween the densities profiles obtained from the different fits 
and the density obtained from the original GPE. The blue 
solid line corresponds to the g-Gaussian approach, while the 
red dashed line corresponds to the traditional Gaussian and 
the green dash-dotted corresponds to the TF approximation. 
As in Fig. IIUI each panel, from top to bottom, corresponds 
to: a) N ^ 10-\ h) N = 10°, c) TV = 10^ d) iV = 10^ and 
e) iV = 10*. 

in the interesting, and more challenging, case of interme- 
diate density N = 10, where the condensed cloud profile 
is neither close to Gaussian nor parabolic, the g-Gaussian 
clearly outperforms both the traditional Gaussian and the 
parabolic Thomas- Fermi approximation by about one or- 
der of magnitude. 

Therefore, the results presented above, are clear evi- 
dence that not only the g-Gaussian is able to provide a 
good variational ODE reduction of the full GPE across 
all the different density regimes, but it is also has a bet- 
ter (or equal) performance than the separate, traditional, 
approximations at either extreme of the density regime. 

VI. CONCLUSIONS 

In this paper we have proposed an efficient variational 
method to investigate the spatio-temporal dynamics of 



magnetically-trapped Bose-condensed gases. To this end 
we have employed a g-Gaussian trial wave-function that 
interpolates between the low- and the high-density limit 
of the ground state of a Bose-condensed gas. Our main 
result consists of reducing the Gross-Pitaevskii equation 
to a set of only three equations: two coupled nonlinear 
ordinary differential equations describing the phase and 
the curvature of the wave-function and a separate alge- 
braic equation yielding the generalized width. On the 
analytical side our equations recover those of the usual 
Gaussian variational approach (in the low-density limit), 
and the hydrodynamic equations that describe the high- 
density case. On the numerical side, we have presented a 
detailed comparison between the numerical results of our 
reduced equations and those of the original GPE for the 
case of periodically driven BEG. We used as indicators of 
the g-Gaussian performance the fitting of the peak den- 
sity, the condensate width, and g(t), a direct comparison 
between the full width at half maximum obtained from 
our variational treatment and that of the original equa- 
tion, and the mean square deviation calculations. All 
these indicators demonstrate that the g-Gaussian vari- 
ational reduction performs very well over all the den- 
sity regimes. Furthermore, we found that a) in the low- 
density regime, the g-Gaussian clearly outperforms the 
traditional Gaussian; b) in the high-density limit the g- 
Gaussian has a similar performance to the Thomas-Fermi 
approximation; c) and that in the intermediate density 
regime the g-Gaussian clearly outperforms both the tra- 
ditional Gaussian and the Thomas-Fermi approximation. 

Finally, let us mention that while our equations are 
one-dimensional they can be easily extended to two- and 
three-dimensional clouds. This avenue is currently under 
investigation and will be reported elsewhere. 
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